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ABSTRACT 

The gap formation induced by a giant planet is important in the evolution of the planet 
and the protoplanetary disc. We examine the gap formation by a planet with a new 
formulation of one-dimensional viscous discs which takes into account the deviation 
from Keplerian disc rotation due to the steep gradient of the surface density. This 
formulation enables us to naturally include the Rayleigh stable condition for the disc 
rotation. It is found that the derivation from Keplerian disc rotation promotes the 
radial angular momentum transfer and makes the gap shallower than in the Keplerian 
case. For deep gaps, this shallowing effect becomes significant due to the Rayleigh 
condition. In our model, we also take into account the propagation of the density waves 
excited by the planet, which widens the range of the angular momentum deposition 
to the disc. The effect of the wave propagation makes the gap wider and shallower 
than the case with instantaneous wave damping. With these shallowing effects, our 
one-dimensional gap model is consistent with the recent hydrodynamic simulations. 

Key words: accretion, accretion discs, protoplanetary discs, planets and satellites: 
formation 


1 INTRODUCTION 


A planet in a protoplanetary disc gravitationally interacts 
with the disc and exerts a torque on it. The torque ex¬ 
erted by the planet dispels the surroundin g gas and forms 
a disc gap along the orbit of the planet llLin fc Papaloizoul 
119791 : iGoIdreich fc Tremainell 19801 '). However, a gas flow into 
the gap is also caused by viscous diffusion and hence the 
gap depth is determined by the balance between the plan¬ 
etary torque and the viscous diffus ion. Accordingly, only a 
large planet can create a d eep ga p _|Lin_&^aMloizo2ll99^; 
Takeuchi et al.lfl996l: I Wardll 19971 : iRafikovI 20021 : Crida etldl 


200611 . 


The gap formation strongly influences the evolution 
of both the planet and the protoplanetary disc in various 
ways. For example, a deep gap prevents disc gas from ac- 
creting onto the plan e t and slows down the planet growth 
dD’Angelo et al.ll2002l : iBate et ai1l2003l : iTanigawa fc Ikomal 
l2007l l. and also changes the pl anetary migration from 
the type I t o the slower type II llLin fc Papaloizoul 1 19861 : 
IWardl 119971 '). Furthermore, a sufficient ly deep gap in- 
hibits gas flow across the gap (lArtvmowicz fc LubowHl996i : 
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iLubow et al.l 1 19991 : iKlevI 119991 1 , which is a possible m echa¬ 
nism for forming an inner hole i n the disc dZhu et al.ll201ll : 
iDodson- Robinson fc Salv^l201ll l. 


Because of their importance, disc gaps induced by plan¬ 
ets have been studied by man y authors, using simp l e one¬ 
dimensional disc mode l s (e.g., TakeuchLeLal, 19961 : IWardl 
Il997l : ICrida et alJ 120061 : ILubow fc D’Ange lo||2006| )~and nu¬ 
merical hydrodynamic simulations (lAr tvmowicz & Lubow 


1994 : lKlevlll999l: 
2013 :lFung et al 


Vanhere et al.|[2004l ; IDuffell fc MacFadven 

I20l4l . One-dimensional disc models pre¬ 


dict an exponential dependence of the gap depth. That is, 
the minimum surface density at the gap bottom is pro¬ 
portional to exp[— A(M P /M,) 2 ], where M p and M, are the 
masses of the planet and the central star, and A is a non- 
dimensional parameter (see also eq. 1381 1. On the other hand, 
recen t hi gh-resolution hydrodynamic simulations done by 
IDuffell fc MacFadvenl (120131 . hereafter DM13) show that the 
gap is much shallower for a massive planet than the pre¬ 
diction of one-dimensional models. According to their re¬ 
sults, the minimum s urface density at the gap is propor- 
tional to (M„/M* 1~ 2 . IVarniere et al.l (l2004lf and lFung et alJ 
J2014ll obtained similar results from their hydrodynamic 
simulations. Its origin ha s not yet bee n clar ified by the one¬ 
dimensional disc model. iFung et al.l (l2014f ) also estimated 
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the gap depth with a “zero-dimensional” analytic model, by 
simply assuming that the planetary gravitational torque is 
produced only at the gap bottom. Their simple model suc¬ 
ceeds in explaining the dependence of the minimum surface 
density of oc (M p /M„)~ 2 . However, the zero-dimensional 
model does not give the radial profile of the surface density 
(or the width of the gap). It is not well understood what kind 
of profile accepts their assumption on the planetary torque. 
Further development of the one-dimensional gap model is 
required in order to clarify both the gap depth and width. 
Such a model enables us to connect the gaps observed in 
protoplanetary discs with the embedded planets. 

One of the problems of the one-dimensional disc model 
is the assumption of the Keplerian rotational speed. The disc 
rotation deviates f rom the K epleria n speed due to a radial 
pressure gradient dAdachi et al.lll976h . When a planet cre¬ 
ates a deep gap, the steep surface density gradient increases 
the deviation of the disc rotation significantly, which affects 
the angular momentum transfer at the gap (see Sections l2.ll 
and !2.2l) . Furthermore, a large deviation of the disc rotation 
can a lso violate the Rayle igh stable condition for rotating 
discs (IChandrasekhaJl939h . A violation of the Rayleigh con¬ 
dition promotes the angular momentum transfer and makes 
the surface density gradient shallower s o that th e Ra yleigh 
cond i tion is only marginal ly satisfied dTanigawa fe Ikomal 
120071: 1 Yang fe Menoul [ 20 T 0 I I. To examine such feedback on 
the surface density gradient, we should naturally include 
the deviation from the Keplerian disc rotation in the one¬ 
dimensional disc model. 

Another simplification is in the wave propagation at 
the disc-planet interaction. The density waves excited by 
planets radially propagate in the disc and the angular mo¬ 
menta of the waves are deposited on the disc by damping. 
This angular momentum deposition is the direct cause of the 
gap formation. Most previous studies simply assume instan- 
tane ous damping of the density wa ves after their excitation 
(e. g.. I Wardll 199 71 : 1 C rid a et alJl2006f) . If the wave propagation 
is taken into account, the angular momentum is deposited 
in a wid er region of the disc, which increa ses the width of 
the gap dTakeuchi et al.lfl996l : lRafikovll2002l l. In a wide gap, 
the disc-planet interaction would be weak because the disc 
gas around the planet decreases over a wide region. Hence, 
we cannot neglect the effect of wave propagation on the gap 
formation. 

I 11 the present paper, we re-examine the gap formation 
by a planet with the one-dimensional disc model, taking 
into account the deviation from Keplerian rotation and the 
effect of wave propagation. To include the deviation from 
Keplerian disc rotation, we modify the basic equations for 
one-dimensional accretion discs, detailed in the next section. 
The effect of the wave propagation is included using a simple 
model. In Section [3] we obtain estimates of gap depths for 
two simple cases. One estimate for a wide gap corresp on ds t o 
the zero-dimensional model proposed bv lFung et al.l (12014) . 
In Sections[4]and[5l we present numerical solutions of the gap 
without and with wave propagation, respectively. We find 
that the gap becomes shallow due to the effects of the devi¬ 
ation from Keplerian rotation, the violation of the Rayleigh 
condition and the wave propagation. With these shallowing 
effects, our results are consistent with the recent hydrody¬ 
namic simulations. In Section [6l we summarize and discuss 
our results. 


2 MODEL AND BASIC EQUATIONS 

We examine an axisymmetric gap in the disc surface den¬ 
sity around a planet by using the one-dimensional model 
of viscous accretion discs. Although the Keplerian angular 
velocity is assumed in most previous studies, we take into 
account a deviation from Keplerian disc rotation in our one¬ 
dimensional model. The deviation cannot be neglected for 
a deep gap, as will shown below. We also assume non-self- 
gravitating and geometrical thin discs. For simplicity, the 
planet is assumed to be in a circular orbit. We also adopt 
simple models for density wave excitation and damping to 
describe the gap formation. 


2.1 Angular velocity of a protoplanetary disc with 
a gap 


The angular velocity, Q, of a gaseous disc around a central 
star with mass M* is determined by the balance of radial 
forces: 


Q. 2 R- 


GM t 

R 2 


1 <9P 2 d 

E OR 


= 0 , 


(1) 


where R is the radial distance from the central star, E is the 
surface density of the disc and P 2 d denotes the vertically 
averaged pressure. On the left-hand side of equation ©, 
the first term represents the centrifugal force on a unit mass 
of the disc. The second and third terms are the gravitational 
force by the central star and the force of the radial pressure 
gradient, respectively. For P 2 d, we adopt the simple equation 
of state P 2 d = c 2 E, where c is the isothermal sound speed. 
Using this equation of state, equation © can be rewritten 
as 


b 2 = (i - 2 v) , 


(2) 


with 


h 2 ( d In E d In c 2 \ 

2 R V dR + dR ) ’ 


( 3 ) 


where 12k = \JGM„ / R 3 is the Keplerian angular velocity, 
and h = c/fl k 

For a disc with no gap, the order of magnitude of 
the n on-dimensional parameter 17 is O {h 2 /R 2 ) (lAdachi et alJ 
EiffiJ), because the term in parentheses in equation m is 
comparable to ~ 1/7?.. On the other hand, if a planet opens 
a deep gap with a width of ~ h, the steep gradient of the sur¬ 
face density increases 17 to 0(h/R). Hence, it also enhances 
the deviation of the disc rotation from the Kepler rotation 
in the deep gap. We neglect the term dc/dR in equation (CTf l 
because the temperature gradient would be small. Neglect¬ 
ing the smaller terms of 0(h/R), we approximately obtain 
dVL/dR as 

dfl _ 3fl K L ft 2 a 2 In El 

dR ~ 27?. [ 3 dR 2 \ ' 

Note that the second term in the parentheses is of order 
unity since d 2 In E/da: 2 ~ 1/h 2 in a deep gap. Therefore, it 
is found that dQ,/dR is significantly altered from the Kep¬ 
lerian value due to the steep gradient of the surface density, 
though the deviation of Q is small (~ h/R). As shown later, 
this deviation promotes radial viscous transfer of angular 
momentum and makes a gap shallower. 
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2.2 Basic equations describing a disc gap around 
a planet 


The equations for conservation of mass and angular momen¬ 
tum are given by 


<9E 1 dF m 

dt + 2ti -R dR 


( 5 ) 


and 


i w + 


1 dFj 
2nR dR 


= jSm+ 2FR Ad ’ 


( 6 ) 


where Fm and Fj are the radial fluxes of mass and angular 
momentum, and j{= R 2 Q) is the specific angular momen¬ 
tum. In equation ©, the source term Sm represents the 
mass accretion rate onto a unit surface area of the disc. The 
accretion of disc gas onto the planet can be included in Sm 
as a negative term. In equation ©, Ad (R) represents the 
deposition rate of the angular momentum from the planet 
on the ring region with radius R. 

To describe the deposition rate Ad, we consider the an¬ 
gular momentum transfer from the planet to the disc. This 
transfer process can be divided into two steps. First the 
planet excites a density wave by the gravitatio nal interac¬ 
tion with the disc (e.g.. IGoldreich fc Tremaine! 1 1980h . Sec¬ 
ond, the density waves are graduall y damped due to the 
disc viscosity or a nonlin ear effect llTakeuchi et al.l Il99fil : 
iGoodman fe Ra fikov 2001). As a result of the wave damp¬ 
ing, the angular momenta of the waves are deposited on the 
disc. If instantaneous wave damping is assumed, the depo¬ 
sition rate Ad is determined only by the wave excitation. In 
Section 12.4.11 we will describe the deposition rate for the 
case with instantaneous wave damping. In Section T2.4.2I we 
will give a simple model of Ad for the case of gradual wave 
damping. 

The ra dial angu lar mom entum flux Fj is given by (e.g., 
iLvnden-Bell fc PringleHl974h 


Fj = jF M - 2 nR 3 vZ^. 


( 7 ) 


The first term is the advection transport by the disc radial 
mass flow, Fm, and the second term represents the viscous 
transport. For the kinetic viscosity, we adopt the a prescrip¬ 
tion, i.e., v = ach (IShakura fc Sunvaevll 19731) . Note that Fj 
does not include the angular momentum transport by the 
density waves in our formulation. 

Equations © © describe the time evolution of the 
three variables E, Fm and Fj with the given mass source 
term Sm, the angular momentum deposition rate from a 
planet Ad and the disc angular velocity fi. Note that Q. de¬ 
pends on dT,/dR, as in equations 0 and 0. 

Next we consider the disc gap in a steady state ( d/dt = 
0). The time scale for the formation of a steady gap is 
approximately equal to the diffusion time within the gap 
width, tdiff = h 2 /v. For a nominal value of a (~ 10 —3 ), 
the diffusion time is roughly given by 10 3 Keplerian peri¬ 
ods, which is shorter than the gr owth time of planets (10 5 7 
yr) ( Kokubo fc Ida 211(11 ). [2002il or the life time of proto¬ 
planetary discs (10 6 7 yr) ( Haisch et al.ll200lh . Hence the 
assumption of a steady gap would be valid. In addition, 
we assume Sm = 0 for simplicity. Although gas accretion 
onto the planet occur s for M p > 10M^ (|Mizunol Il98fll : 
iKanagawa fe Fuiimotd l2013li , the assumption of Sm = 0 


would be valid if the accretion rate onto the planet is smaller 
than the radial disc accretion rate, Fm. 

Under these assumptions, equation © shows that Fm 
is constant. Equation © yields 

/> OO 

Fj = Fj(oo) - / Ad dR', (8) 

Jr 

where Fj(oo) is the angular momentum flux without the 
planet. From equations © and ©, we obtain 


HD f°° 

jF M - 2 -kR 3 T,v— = Fj(oo) - / A d dR 1 . (9) 

an J R 

Equation (© with a constant mass accretion rate describes 
a steady disc gap around a planet for a given Ad. Since 
dti/dR is given by equation 0, equation 0 is the second- 
order differential-integral equation. Note that equation © 
is derived from equation © and indicates the angular mo¬ 
mentum conservation. 

By differentiating equation ©, we obtain a rather fa¬ 
miliar expression for the mass flux: 


Fm 



_d_ 

dR 




+ Ad 


( 10 ) 


Note that this expression is valid only for the steady state. 
In a time-dependent case, equation m should include the 
term — 2nRS(dj/dt) in the parentheses. As a boundary con¬ 
dition, the disc surface density should approach its unper¬ 
turbed values at both sides of the gap far from the planet. 

Here, we also consider the unperturbed surface density. 
In the unperturbed state, S2 can be replaced by S7k, by ne¬ 
glecting the smaller term O [h 2 / R 2 ) (see equations [5] and 
ED- Furthermore, setting Ad = 0 in equation (|9|l . we obtain 
the unperturbed surface density Eo as 

3nR 2 un K T, 0 (R) = -F 2 B K F M + Fj(oo). (11) 


Thus Eo is given by 

v _ fi _ -^(00) \ 
3ni/ y R 2 QkFm ) 


( 12 ) 


This agrees with the well -kn own solution for stea dy viscous 
accretion discs ('e.g.. ILvnden-Bell fc Pringlelll974h . 


2.3 Rayleigh condition 


For a deep gap around a large planet, the derivative of the 
angular velocity deviates significantly from the Keplerian 
velocity, as shown in Section 12.11 A sufficiently large de¬ 
viation in fi violates t he so-called Rayleigh stable condi¬ 
tion of dj/dR ^ 0 (see. IChandrasekharill96If) . Such a steep 
gap is dynamically unstable, which would cause a strong 
angular momentum transfer, lessening the steepness of the 
gap. This would make the unstable region marginally stable 
(i ,e.,dj/dR = 0). 

Using equation ©. we give dj/dR. as 


dj_ 

dR 


— — 7? D fl 


Kp 


1 + h v 


d 2 In E \ 
~dR 2 ~) ’ 


(13) 


where the suffix p indicates the value at R = R p - this suffix 
is also used for other quantities. Hence, using the second- 
derivative of the surface density, the m arginally stable con- 
dition dj/dR = 0 can be rewritten as dTanigawa fe Ikomal 
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|2007j). 


, 2 d 2 In E 
h _ 

p dR 2 


- 1 . 


(14) 


Actually, around a sufficiently large planet, equation © 
gives hpd 2 In T,/dR 2 < —1 in some radial regions. In such 
unstable regions, we have to use equation m instead of 
equation ©• 

The breakdown of equation m indicates that the flux 
Fj of equation 0 cannot transport all of the angular mo¬ 
mentum deposited by the planet. In a real system, however, 
the instability would enhance the angular momentum flux, 
which keeps the gap marginally stable. The enhancement of 
Fj can be considered to be due to an effective viscosity u e g 
enhanced by the instability. Since such an effective viscosity 
restores equation m , r'efi in the unstable region is given by 


Z'eff = 


—jFm + Fj(oo) — A d dR' 

47rF 2 EfI 


(15) 


where we use the relation dQ/dR = —2 fl/R obtained from 
the marginally stable condition. Furthermore, by using u e g 
instead of v, equation 0 gives the enhanced angular mo¬ 
mentum flux in the unstable region. 

The Rossby w ave instability may be important for the 
gap f ormation (e.g.. lRichard et al.ll2013l : IZhu et al.ll2013l : iLinl 
1201411 . As well as the Rayleigh condit ion, the Rossb y wave in¬ 
stability relates to the disc rotation dLi et al ■2Q03). Because 
it can occur before the Rayleigh condition is violated, how¬ 
ever, the Rossby wave instability may suppress the surface 
density gradient more than the Rayleigh condition. For sim¬ 
plicity, we include only the Rayleigh condition in the present 
study. A further detail treatment including the Rossby wave 
instability should be done in future works. 


2.4 Angular momentum deposition from a planet 

In the disc-planet interaction, a planet excites density waves 
and the angular momenta of the waves are deposited on 
the disc through their damping. The angular momentum 
deposition rate Ad is determined by the later process. Firstly, 
we will consider the deposition rate Ad in the case with 
instantaneous wave damping. In this case, the deposition 
rate is governed only by the wave excitation. Next, taking 
into account the wave propagation before damping, we will 
model the deposition rate in a simple form. 


2.4-1 Case with instantaneous wave damping 

Under the assumption of instantaneous wave damping, the 
angular momentum deposition rate Ad(-R) is equal to the 
excitation torque density A ex (.R), which is the rate at which 
a planet adds angular momenta to density waves per unit 
radial distance at R. That is, 

A d = A ex . (16) 

At a position far from the planet, the ex citation torq ue den¬ 
sity is given by the WKB formula (e.g.. I Wardll 1986h as 

AS™ -±c*ni e (0 1 («^) 4 . («) 


where C = (2 5 /3 4 )[2A 0 (2/3) + A'i(2/3)] 2 /tt ~ 0.798 and 
Ki denote the modified Bessel functions. The sign of equa¬ 
tion (1 1711 is positive for R > R p or negative for R ^ R p . In 
the close vicinity of the planet, |A — R p \ < h p , on the other 
hand, the WKB formula is overestimated. Thus, we model 
the excitation torque density A ex with a simple cutoff as 

A (A-- for \R — A p | > h p A, 

\0 for \R — A p | < h p A. v ' 

The cut-off length h p A is determined so that the one¬ 
sided torque T (= A ex dR) agrees with the result of the 

linear theory for realistic discs (iTakeuc hi fc Mivamal 1 19981 : 

iTanaka et al.ll2002l : iMuto fe Inutsukall2009l l. Then we obtain 

A = 1.3. 

Note t hat the WK B formula is derived for discs with no 
gap. IPetrovich fe Rafikovl (I2012T ) reported that the torque 
density is altered by the steep gradient of the surface den¬ 
sity because of the shift of the Lindblad resonances. For 
simplicity, however, we ignore this effect in the present pa¬ 
per. Hence, in our model, the excitation torque density 
A ex is simply proportional to the disc surface density at 
R, E (R), and is independent of the surface density gra¬ 
dient even for deep gaps. For a large planet with a mass 
of Mp/M* > ( h p /R p ) 3 ; furthermore, the non-li near effect 
would not be negl igible for wave excitation dWardl Il997l : 
IMivoshi et al.ll 19991) . This non-linear effect is also neglected 
in our simple model. 


2.4-2 Case with wave propagation 


When wave propagation is included, the angular momen¬ 
tum deposition occurs at a different site from the wave 
excitation and equation m is not valid. In this case, 
the angular momentum deposition is also governed by the 
damping of the waves. Although the wave damping has 


been examin ed in previous studies (e.g., ITakeuchi et al 


199C; IKorvcanskv fe Papaloizoulll996l : [Goodman fe Rafikov 


200111 . it is not clear yet how the density waves are damped 
in a disc with deep gaps. In the present study, therefore, 
we adopt a simple model of angular momentum deposition, 
described below. 

Since the waves are eventually damped in the disc, the 
one-sided torque (i.e., the total angular momentum of the 
waves excited at the outer disc in unit time) is equal to the 
total deposition rate in the steady state. That is, 


P OO POO 

T= AexdA' = / A d dR’. 

•J Rrt J Rn 


(19) 


Using the one-sided torque, the angular momentum deposi¬ 
tion rate can be expressed by 


A d = ±Tf(R), (20) 

where the distribution function f(R) satisfies f(R)dR = 
1, and the sign is the same as in equation m- As a simple 
model, we assume a distribution function f(R) given by 

for x d h p - ^ < \R - A p | < x d h p + ^, 
otherwise. 



( 21 ) 
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In this simple model, the non-dimensional parameter xa de¬ 
termines the position of the angular momentum deposition 
and the parameter Wd represents the radial width of the de¬ 
position site. The waves propagate from the excitation site 
to the deposition site around |x| = xa■ Since the density 
waves propagate away from the planet, the deposition site 
is farther from the planet than the excitation site. The pa¬ 
rameter xa should be consistent with this condition. 

In the case with wave propagation, we use equa¬ 
tions (|20l> and m to obtain the gap structure with equa¬ 
tion ©• It should be noted that T in equation EQ| depends 
on the surface density distribution through the definition 
of equation G9, because A ex is proportional to E. These 
coupled equations are solved as follows. First, we obtain the 
surface density distribution with equation © for a given T. 
Next, we determine the corresponding mass of the planet 
from equation m, using the obtained surface density. 


2.5 Local approximation and non-dimensional 
equations 


The typical width of a disc gap is comparable to the disc 
scale height and much smaller than the orbital radius of 
the planet. Thus it is convenient to use the local coordinate 
defined by 


x = 


R — R p 

tip 


( 22 ) 


Note that the suffix p indicates the value at R = R p . 

We adopt a local approximation in which terms propor¬ 
tional to hp/Rp and higher order terms are neglected. From 
equations © and ©, the deviation in fl from fix is given 
by 


/ip^Kp c?lnS 

n ~^ = 2Rp dx 


(23) 


Rp given by equation m■ Dividing equation ESI by 
37rf?.pr' p Eo(f?p)flKp and using equation (1111) . we obtain a 
non-dimensional form: 


1 d 2 In s \ 
3 dx 2 J 


s 


1 , 
= 1 - - J A d dx', 


(27) 


where Ad is the non-dimensional angular momentum depo¬ 
sition rate defined by 


Ad 


Ad ftp 

nRpU p T,o(Rp)tiKp ' 


(28) 


The marginally stable condition can be rewritten as 

d 2 In s _ 
dx 2 


(29) 


This equation is used instead of equation m in the 
Rayleigh unstable region. 

The non-dimensional excitation torque density, A ex , is 
defined by 


A ex fo p __ I ±K—s(x) for |x| > A, 

7rRpt'pEo(Rp)nK P |q f or ^ a. 

(30) 


where the non-dimensional parameter I\ is given by 



(31) 


In the above, we use u p = a/ipflxp. In our model, the pa¬ 
rameter K is the only parameter that determines the gap 
structure for the instantaneous damping case. 

In the case with instantaneous wave damping, the angu¬ 
lar momentum deposition rate is given by Ad = A ex feci. 1161 1. 
In the case with wave propagation, equation ©3 gives 


A d (x) = fhpf{R v + hpx), (32) 


and is proportional to h p /R p . Thus, the disc angular veloc¬ 
ity ft is replaced by the angular velocity of the planet flxp 
under the local approximation, and the specific angular mo¬ 
mentum j also is given by R 2 Qk p - As for the derivative 
dQ/dR, we cannot neglect the deviation from the Keplerian 
value. Equation © yields 


dQ 3Dk p f ld 2 lnE\ 
dR = 2 Rf \ ~ 3 dx 2 J ' 


(24) 


Equation © can be rewritten in the local approxima¬ 
tion as 


R 2 HkpFm+ 37ri?p Up EDk p 


/ 1 d 2 In E \ 

v 3~chp- J 


POO 

= Tj(oo) — / Kahpdx . 
J X 


(25) 


Because of the local approximation, equation (ESI cannot 
be applied for the wide gap formation. If the half width of 
gap is narrower than about l/3i? p , equation (1251) would be 
valid. 

Here we introduce the non-dimensional surface density, 
s, defined by 


E 

So(-Rp) ’ 


(26) 


where Eo(f? p ) is the unperturbed surface density at R = 


where the non-dimensional one-sided torque, T, is given by 

POO 

T — K —^s(x)dx. (33) 

J a r, 

Note that the deposition rate Ad includes two parameters Xd 
and Wd, in addition to K. In this case, we solve equation m 
for a given value of T. Then, we can obtain K by substituting 
the solution s(x) into equation (1331) . as mentioned at the end 
of the last subsection. In order to obtain the solution for a 
certain K, we need an iteration of the above procedure with 
trial values of T. 

The boundary conditions of equations (1271) and (1291) are 

s = 1, at x = ±oo. (34) 

Under the local approximation, the surface density has a 
symmetry of s(x) = s(—x), since both the above basic equa¬ 
tions and the deposition rate are symmetric. 


3 ESTIMATES OF GAP DEPTHS FOR 
SIMPLE SITUATIONS 

3.1 Case of the Keplerian discs 

Before deriving the gap solution in our model described in 
Section © we examine the gaps for two simple situations. 
First, we consider a disc with Keplerian rotation, as assumed 
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in previous studies. Neglecting the deviation in dfl/dR from 
the Keplerian (i.e., the term of d 2 In s/dx 2 ) in equation (1271) . 
we have 

1 r°° 

s = 1 — - / Ad dx' (35) 

3 J X 

Here we also assume instantaneous wave damping and adopt 
Aa = A ex (eq. [30] )• Differentiating equation (1351) . we obtain 


dins 

dx 


A——y A 

for 

3* 4 

0 

for 


M > A, 
1*1 A, 


(36) 


Hence, we obtain the surface density in the Keplerian discs 
with the instantaneous wave damping as 


s(x) = 


exp | 


for 

exp | 

:-&*) 

for 


for 1*1 > A, 


(37) 


Using equation m, C = 0.798 and A = 1.3, the minimum 
surface density, s m in, is 


Smin — GXp 


—0.040a 




5 / JWp X 2 

M* 


(38) 


This solution is almost the sa me as that in the previou s 
one-dimensional gap model (e.g., Lubow & D’A ngelo 120061 ). 

For a very large K, the Rayleigh c ondition is vi olated 
and equations (I3Z1) and (EH) are invalid. iTanigawa fe Ikomal 
(l2007ll obtained the gap structure in Keplerian discs, in¬ 
cluding the Rayleigh condition. Their solution is described 
in Appendix m In Appendix m we also derive gap solutions 
in Keplerian discs, taking into account the wave propagation 
with the simple model of equation (1201) and ED. 


3.2 Case of the wide-limit gap 

Next, we conside r a situation implied by the zero-dimension 
analysis done bv lFung et al.l (( 2014h . which assumes that the 
wave excitation occurs only a the gap bottom. This assump¬ 
tion would be valid if the gap bottom region is wide enough. 
Hence, we call this situation ’wide-limit gap’ case. Since the 
density waves are excited at the gap bottom with s ~ s m i n , 
the one-sided torque of equation (1331) is simply given by 

T = ^3 ATs min ~ 0.121/Cw (39) 

Using equation (1271) . we can estimate s m i n of the wide-limit 
gap. The right-hand side of equation ED can be rewritten 
as 1 — T/3 at x = 0. In the left-hand side of equation (1271) . 
moreover, we can neglect the term d 2 In s/dx 2 when a flat- 
bottom gap is assumed. Then, the relation between s m i n and 
T is obtained as 

Smin = 1 - J. (40) 

Equations (1391) and ED yield 

Smin= i + o.o4o/r (41) 

For a large K, s m i n given by equation ED is proportional to 
1/A'. This result agrees with the zero-dimensional model by 


iFung et alj (120141 (1*1. In the zero-dimensional model, the min¬ 
imum surface density is estimated from a balance between 
the planetary torque and the viscous angular momentum 
flux outside the gap. Such a balance is also seen from equa¬ 
tion (1401) (and eq. |27[). The first and second terms in the 
right-hand side of equation (BED correspond to the viscous 
angular momentum flux outside the gap and the planetary 
torque and the left-hand side is negligibly small for a large 
K. 

With their hydrodynamic simulations for K < 10 4 , 
DM13 derived a similar resuh@, 


29 1 

29 + K ~ 1 + 0.034A ' 


(42) 


It is found that equations OD and ED are consistent with 
each other. Note that these minimum surface densities are 
much larger than that of the Keplerian disc (eq. [38]) for a 
large K because equation (1391) is not accepted in the Ke¬ 
plerian solution. The wide-limit gaps assume that all the 
waves are excited in the bottom region with s ~ s m i n , i.e., 
equation ED- In Sections [4] and [5] we will check whether 
or not this assumption is valid, by comparing it with our 
one-dimensional solutions. 


4 GAP STRUCTURE IN THE CASE WITH 
INSTANTANEOUS WAVE DAMPING 


4.1 Linear solutions for shallow gaps 

Here we present the numerical solution of the gap in the case 
with instantaneous wave damping (i.e., Aa = A ex ). 

First, we consider the case with a small K in equa¬ 
tion ED^ in which Aa is proportional to K. This case corre¬ 
sponds to a shallow gap around a small planet. Since |s — 1| 
is small, it is useful to express the solution as 


s = exp(Ay), 


(43) 


or s = 1 + Ky. As seen in the next subsection, the former 
expression is better for an intermediate K (~ 10). Substi¬ 
tuting equation ED into equation ED with equation ED, 
we can expand it into a power series of Ky. The first order 
terms give the linear equation of y: 


d*y 
dx 2 


- 3 y 


-F 


C 

3|*| 3 

C 

3A 3 


for |*| > A, 
otherwise, 


(44) 


where the sign in the right-hand side is negative for * > 0 
and positive for * ^ 0. Equation ED is an inhomogeneous 
linear differential equation, and can be integrated with the 
boundary conditions of equation ED- We do not need to 
take care of the Rayleigh condition in the shallow gaps. A 
detailed derivation of the linear solution is described in Ap¬ 
pendix [C] 

Fig.[Dt shows y, which can be converted into the surface 
density s by equation ED- In these shallow gaps, the gap 
depth is almost the same as for the Keplerian case, though 
our model gives a smooth surface density distribution. 


1 In the notation of lFung et alj d2014T ). K is given by q 2 /(a[h/ r] 5 ) 

2 In the notation of DM13, K is given by Af _1 (M s h/M p ) 2 » _1 . 
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Figure 1. Linear solution for y (a) and dy/dx (b). The surface 
density and angular velocity are given by y and dy/dx with equa¬ 
tions 03 and iTTsT i . respectively. The dashed line is the solution 
for the Keplerian disc. 


-4-3-2-101234 

X 

Figure 2. Surface density distributions for K = 50 (a) and 200 
(b). The red line is the exact solution (see text). The chain line 
is the linear solution given by equation S3 and the dashed line 
is the solution for the Keplerian case (eq. E2). The dotted line 
represents the minimum surface densities for the wide-limit gap 
given by equation tni 


Fig.[T|) shows the derivative of y which is related to AS2, 
as 


An = n-n K =K^^f-, (45) 

ZHp ax 

using equations (1231) and (1431) . The absolute value of Atl 
attains a maximum at |x| ~ 1.5. The second-order derivative 
of y gives the shear, dQ/dx , as 


dfl d^K F K d 2 y\ 

dx dx \ 3 dx 2 ) ’ 


(46) 


as seen from equation (1241) . At |x| > 1.5, the shear mo¬ 
tion is enhanced compared to the Keplerian case, because 
d 2 y/dx 2 < 0. Since the shear motion causes viscous angu¬ 
lar momentum transfer, this enhancement makes the sur¬ 
face density gradient less steep compared with the Keplerian 
case, as shown in Fig. [T}i. 


4.2 Nonlinear solutions for deep gaps 

Next we consider deep gaps around relatively large plan¬ 
ets. In this case, we numerically solve the non-linear equa¬ 
tion m with the Rayleigh condition. We call the obtained 
non-linear solution the “exact” solution. 


At regions far from the planet, the surface density per¬ 
turbation is rather small and the linear approximation is 
valid. Thus, we adopt a linear solution at |x| > 10. Note that 
this linear solution has different coefficients for the homoge¬ 
neous terms from those in Section|4Tjsee Appendix 0. The 
coefficients of the homogeneous solution are given to satisfy 
the boundary conditions of equation (1341) . At |cc| ^ 10, we 
integrate equation (|27l) with the fourth-order Runge-Kutta 
integrator. In the Rayleigh unstable region, the surface den¬ 
sity is governed by the marginally stable condition (eq. 1291 1. 

instead of equation 123- 

Fig. [2]shows the surface density distributions of the ex¬ 
act solutions for K = 50 (a) and 200 (b). If we assume a disc 
with h p /R p = 0.05 and a = 10 —3 , these cases correspond to 
M p = l/8Mj and l/4A/j respectively, where Mj is the mass 
of Jupiter. For comparisons, the Keplerian solution (eq. 1371 1 
and the linear solution with equation (03 are also plotted. 
For K = 50, the linear solution almost agrees with the exact 
solution, while it is much deeper than the exact solution for 
K = 200. For I\ = 200, the Keplerian solution has a much 
smaller s m i n than the exact solution. 

Fig. [3 illustrates the angular velocities (a) and specific 
angular momenta (b) for the exact solutions for K = 50 and 
200. Similar to the linear solution in Fig.Q] the shear motion 
is enhanced at \x\ > 1.4. This enhancement of the shear mo- 
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Figure 3. (a) Deviation from Keplerian disc rotation and (b) 

specific angular momentum, for K = 50 (dashed) and 200 (solid). 
The filled circles indicate the edge of the marginally stable region 
for the Rayleigh condition. 



Figure 4. Shear of exact solutions for K = 50 (dashed) and 200 
(solid). 


tion is also seen in Fig. [4] The enhancement promotes the 
angular momentum transfer and makes the surface density 
gradient less steep. For K = 200, the Rayleigh condition is 
violated. In the unstable region, the marginally stable condi¬ 
tion further reduces the surface density gradient. This makes 
the gap much shallower than for the Keplerian solution, as 
seen in Fig. [2jn 

We also plot the minimum surface densities, s m i n , of the 
wide-limit gap (eq. [IT]) in Fig. [2] The wide-limit gap gives 
a much larger s m in than the exact solution for I\ = 200. 



X 

Figure 5. Excitation torque density given by equation (1301) for 
K = 200. The two vertical lines indicate the positions with s = 
10s • 

±uo min • 


In the wide-limit gap, it is assumed that the density waves 
are excited only at the gap bottom with s ~ s m in- Fig. [5] 
shows the excitation torque density given by equation m 
for the exact solution with K = 200. This torque density in¬ 
dicates that the waves are excited mainly in the region with 
s > 10s m i n . Thus the assumption of wave excitation at the 
gap bottom is not valid in this case. Since wave excitation 
with a larger s increases the one-sided torque, this can ex¬ 
plain why the gap of the exact solution is much deeper than 
the wide-limit gap in Fig. [5] Note that this result for the 
wave excitation is obtained in the case of instantaneous wave 
damping. The effect of the wave propagation can change the 
gap width and the mode of wave excitation, as seen in the 
next section. 


4.3 Effect of the Rayleigh condition 


We further examine the effect of the Rayleigh condition on 
the gap structure. Fig. [6] shows the surface densities (a) 
and specific angular momenta (b) for the exact solution and 
the solution without the Rayleigh condition. The solution 
without the Rayleigh condition has unstable regions with 
dj/dx < 0 (i.e., 1.4 < \x\ < 3.1). This comparison between 
these two solutions directly shows how the Rayleigh con¬ 
dition changes the gap structure. The Rayleigh condition 
increases s m i n by a factor 6 for K = 200. This is because the 
marginal condition of d 2 In s/dx 2 ^ — 1 keeps the surface 
density gradient less steep and makes the gap shallow. 

It can be considered that the marginally stable state is 
maintained by v e s of equation (USD- The non-dimensional 
form of equation m is given by 


v e g _ 3 - A d dx' 

v 4s 


(47) 


Fig. |T] shows u e ff in the unstable region for K = 200. The 
effective viscosity is twice as large as the original value at 
x = 1.8. This enhancement of the effective viscosity causes 
the shallowing effect in Fig. [SJi. 

In Fig. EK we also plot the surfa ce density distribu¬ 
tion given by Tanigawa fc IkorTial J2007I) (hereafter TI07), in 
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Figure 6. (a) Surface density distribution and (b) specific an¬ 
gular momentum distribution, for K = 200. The red line indi¬ 
cates the exact solution. The green line is the solution with¬ 
out the Rayleigh condition (see text). The chain line in (a) 
denotes the surface dens ity distribution given by the model of 
iTanigawa Sz Ikomal (|2007t) (eq. [A3] . TI07). 



Figure 7. Effective viscosity z/ e ff of the exact solution with K = 
200 . 

which the Rayleigh condition is taken into account (for de¬ 
tails, see Appendix Their model gives a shallower gap 
than our exact solution. This is because a very steep sur¬ 
face density gradient in the Keplerian solution is suppressed 
by the Rayleigh condition to a greater extent than in our 
model. 



Figure 8. Minimum surface densities, s m i n , for the exact solu¬ 
tion (red line) and the solution without the Rayleigh condition 
(green line). The dashed line is s m ; n in the Keplerian case. The 
chain, dotted and solid lines denote .s ln j n given by the model of 
TI07, the wide-limit gap (eq. 1411 1 and the empirical relation of 
DM13 (eq. [42]) respectively. 

We also show that the Keplerian solution by TI07 does 
not satisfy the angular momentum conservation. The Kep¬ 
lerian solution without the Rayleigh condition (eq. [33) is 
derived just from equation (1351) (or eq. EH), which is orig¬ 
inated from equation ©. In this solution, thus, the angu¬ 
lar momentum conservation is satisfied. However, when the 
Rayleigh condition is violated, the marginal stable condition 
(eq. [29]) is used instead of equation (1371) . Because of this, 
the surface density at the flat bottom of TI07’s solution does 
not satisfy equation (1351) or the angular momentum conser¬ 
vation, either. This violation is resolved in our formulation 
because our exact solution always satisfies equation OB out¬ 
side the Rayleigh unstable region 0 

4.4 Gap depth 

Fig- [Hlshows the minimum surface densities, s m in, as a func¬ 
tion of K for the exact solutions. For comparison, we also 
plot s m in for the solutions without the Rayleigh condition 
and the Keplerian solutions. These solutions give deeper 
gaps than the exact solution, similar to the result of Sec¬ 
tion m . It is found that the shallowing effect due to the 
Rayleigh condition becomes significant with an increase in 
K. This is because the Rayleigh condition is violated more 
strongly for large K. 

In Fig. [HJ on the other hand, the exact solution is much 
deeper than DM13’s results and the wide-limit gap, though 
the latter two cases agree well with each other. The model of 
TI07 also gives much deeper gaps than DM13. These com¬ 
parisons indicate that in the case with instantaneous wave 
damping, our exact solution cannot reproduce the hydrody¬ 
namic simulations of DM13. This difference in the gap depth 

3 By introducing the effective viscosity of equation m and mul¬ 
tiplying the LHS of equation 1271) by equation (1271) is 

recovered in the Rayleigh unstable region. 
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from DM13 is likely to be due to the fact that the assump¬ 
tion of the wide-limit gap is not satisfied in the case with 
instantaneous wave damping (see Fig. El). In the next sec¬ 
tion, we will see that the effect of wave propagation widens 
the gap and makes the assumption of the wide-limit gap 
valid. 


5 EFFECT OF DENSITY WAVE 
PROPAGATION 

In this section, we consider the effect of wave propagation. 
Wave propagation changes the radial distribution of the an¬ 
gular momentum deposition. A simple model of angular mo¬ 
mentum deposition rate altered by wave propagation is de¬ 
scribed in Section 12.4.21 Using this simple model, we solve 
equation CD with the Rayleigh condition in the similar way 
to the previous section. At the region far from the planet 
(i.e., |*| > 10), we use the linear solution to equation (1C 1 1) 
with g(x) = 0 in this case. 

5.1 Gap structure for K = 200 

Fig. E] illustrates the surface densities (a) and the excita¬ 
tion torque densities (b) of the exact solutions in the case 
with wave propagation. The angular momenta of the excited 
waves are deposited around |*| = Xd in our model. A large 
Xd indicates a long propagation length between the excita¬ 
tion and the damping. The parameter K is set to 200. For 
an increasing x d , the gap becomes wider and shallower. The 
gap width is directly governed by the position of the angu¬ 
lar momentum deposition. For xa = 3 and 4, the gap depths 
are consistent with the wide-limit gap (and also DM13). For 
Xd = 4, the density waves are excited mainly at the bot¬ 
tom region with s ~ s m i n , as seen in Fig. [9};>. Moreover, 
for Xd = 3, a major part of the wave excitation occurs at 
the bottom. That is, the assumption of the wide-limit gap is 
almost satisfied for the solutions with Xd = 3 and 4. This ex¬ 
plains why the gap depths are consistent with the wide-limit 
gap for these large sa¬ 
lt is also valuable to compare the gap width with hy¬ 
drodynamic simulations. DM13 performed a simulation for 
the case of M p = l/4A/j (2M s h in their notation), a = 10 -3 
and hp/R p = 0.05. This case corresponds to K = 200. In 
this simulation, they found that the gap width is about 6 h v , 
assuming that these gap edges are located at the position 
with E = (l/3)Eo(Rp) (i.e., s = 1/3). If we adopt the same 
definition of the gap edge, the gap widths of our exact solu¬ 
tions with Xd = 3 and 4 are 6.1h p and 7.7 h p , respectively. 
Hence, if we take into account the wave propagation and 
adopt Xd = 3-4, our exact solution can almost reproduce 
both of the gap width and depth of the hydrodynamic sim¬ 
ulations by DM13, for K = 200. 

It should be also noted that, for Xd = 2, the wave excita¬ 
tion mainly occurs at |x| > x d (80% of the excitation torques 
come from this region). However, the deposition site should 
be farther from the planet than the excitation site because 
the density waves propagate away from the planet. Thus, the 
case with Xd = 2 does not represent a realistic wave prop¬ 
agation. From now on, we judge that our simple model for 
the wave propagation is valid if more than half of the one¬ 
sided torque arises from the excitation at |x| < Xd■ In the 
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Figure 9. Surface densities (a) and excitation torque densities 
(b) in the case with the wave propagation for K = 200. The 
green, blue and red lines denote the solutions with x& = 2, 3 and 
4 respectively. The parameter is set to h p . The gray dashed 
line is the surface density in the case with instantaneous wave 
damping. The dotted line in (a) represents the minimum surface 
density for the wide-limit gap given by equation HD , i.e. s min = 
0.109. 


case with Xd = 3 or 4, the excitation at |x| < Xd contributes 
55% or 78% of the one-sided torque, respectively. 

In Fig. 1101 we check the effect of the width of the depo¬ 
sition site, Wd, for * d = 3 and K = 200. It is found that the 
width Wd has only a small influence on the gap structure. 

We show that the deviation from the Keplerian rota¬ 
tion is also important in the case with wave propagation. 
In Fig. 1111 we plot the solution with the Keplerian rotation 
and our exact solution. The Keplerian solution is derived 
from equation (1351) with the angular momentum deposition 
model (eqs. [20] and 1211 1. When the Rayleigh condition is 
violated, the marginal stable condition (eq. [29]) is used. A 
detail derivation of this solution is described in Appendix iBl 
In the Keplerian solution of Fig. 1111 the Rayleigh condition 
is violated over the whole region of the angular momentum 
deposition. Then the minimum surface density is given by 
equation (IB3I) . which is much larger than our solution and 
equation OD. Because equation (IB3I) does not satisfy equa¬ 
tion (13511 , the Keplerian solution does not satisfy the angular 
momentum conservation, as pointed out in Section 14.31 On 
the other hand, in the zero-dimension analysis bv lFung et al.l 
d2014ll (or in eq. ED), Smin is estimated from a balance be- 
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Figure 10. Gap structures for w& = 0.1 h p (green), h p (blue) and 
1.5/ip (light-blue). The parameters K and x& are set to 200 and 
3, respectively. The gray dashed line indicates the solution in the 
instantaneous damping case and the dotted line is the minimum 
surface density of the wide-limit gap (eq. 1411 ). 



x 

Figure 11. The Keplerian solution in the case of wave prop¬ 
agation for K = 200, x ( [ = 4 and w ( i = h v (bule line). For 
comparison, the exact solution (red line) is also plotted. 


tween the planetary torque and the viscous angular momen¬ 
tum flux (i.e., from the angular momentum conservation). 
Because of this difference, the Keplerian solution gives a 
much shallower gap than the estimation in equation 03. 
Note that because our exact solutions are given by equa¬ 
tion 03, the balance between the planetary torque and the 
viscous angular momentum flux is always satisfied in our 
solutions. Hence, our soltuions always satisfy the angular 
momentum conservation and gives a similar s m i n to equa¬ 
tion GD for a sufficiently wide gap. 


5.2 Dependences of the gap depth and width on I\ 

Fig- HU shows the minimum surface densities s m in as a func¬ 
tion of the parameter K similar to Fig. [8] but the effect of 
the wave propagation is included in this figure. In this figure, 
we also show the dependence on the parameter Xd while Wd 
is fixed at h p since Wd does not change the surface density 


Mp/Mj (« = 10 3 , hp/Rp = 0.05) 

0.02 0.05 0.1 0.2 0.5 1 



Figure 12. Minimum surface densities, s m in, of the exact solu¬ 
tions in the case with the wave propagation for Xd = 3 (blue), 4 
(red) and 6 (green). The parameter Wd is h p . We also plot results 
by DM13 (eq. 1421 . solid line) and the wide-limit gap (eq. 1111 . 
dotted line) and the exact solution with instantaneous damping 
(gray line). The dashed lines indicate exact solutions with unre¬ 
alistic wave propagation. 


distribution much (see Fig. II 1)1) . At K = 200, as also seen 
in Fig. [9] our exact solutions reproduce the gap depth of 
DM13 (or the wide-limit gap) for Xd ^ 3. At K = 1000, on 
the other hand, a larger Xd (^ 6) is required for agreement 
with DM13. That is, with an increase of K, a large Xd is 
necessary for the values of the depth of the hydrodynamic 
simulations to agree. Note that the dashed lines in Fig. Q3] 
represent the cases of unrealistic wave propagation, in which 
more than half of the one-sided torque is due to the excita¬ 
tion at |*| > Xd, as for the case of Xd = 2 in Fig. [3 At large 
K, a large Xd is also required for realistic wave propagation. 

We also show the Keplerian solution with the wave 
propagation of Xd = 6 (see Appendix P. For K > 30, •Smin 
is given by equation (lB3l) and independent of K because of 
the Rayleigh condition, as seen in Fig. [9] This unrealistic 
result in the Keplerian solution is related with the violation 
of the angular momentum conservation, as pointed out in 
Section S3] (and see also Appendix P- 

To check which Xd is preferable, we also compare the 
gap width with the hydrodynamic simulations. Fig. ll3l shows 
the gap width of the exact solutions as a function of K. 
Similar to DM13, the gap edge is defined by the position 
with s = 1/3. In this definition, the gap width is roughly 
given by twice Xdh for our exact solutions with K > 50. 
Note that this definition is useless for K < 50 because of 
shallow gaps with s m i n > 1/3. The dashed lines represent the 
cases of unrealistic wave propagati on, sim ilar to Fig. 1121 The 
results of DM13 and Varniere et al.1 (120041 ) are also plotted in 
Fig. 1131 IVarniere et alJ ( 20041 ) also performed hydrodynamic 
simulations of gap formation for M p /M* = 10 4 - 2 x 10 3 , a = 
6x 10 -2 - 6x10 -5 and h p /R p = 0.04 (i.e., K = 600-6x10 s ). 
Their gap depths almost agree with DM13’s relation. For 
K < 300, our exact solu tions wi t h Xd = 3 and 4 agree with 
the results of DM13 and I Varniere et ail (|2004l ). respectively. 
For K > 300, on the other hand, the widths obtained by 
IVarniere et ali (120041 ) are wider than those given by DM13. 
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M p /Mj (a = 10 -3 , h p /R p = 0.05) 
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Figure 13. Gap widths of our solutions for = 3 (blue), 4 (red) 
and 6 (green). The edge of the gap is defined by the position with 
s = 1/3, in the same way as defined by DM13. The parameter Wd 
is set to hp. The dashed lines represent solutions with unrealistic 
wave propagation, similar to Fig. 1121 The circles indicate the gap 
widths obtained by DM13 (Fig. 6 of their p aper ), and the triangles 
show the gap widths by IVarniere et al.l lIlOOlI l (twice Aiooo in 
table 1 of lVarniere et all) . 


Ou r exact solution w ith Xd = 6 agrees with the widths of 
IVarniere et al.l (12004 ). while widths of DM13 correspond to 
our solutions of unrealistic wave propagation. The preferable 
Xd cannot be determined only by this comparison, and we 
still have a large uncertainty in the preferable value of Xd- 

The diffe rence of widths between DM13 and 
IVarniere et al.l (l2004l l would be caused by different param¬ 
eters in their simulations (e.g., the disc viscosity, spatial 
resolution and width of a computation domain). However, 
the origin of the difference is sti ll uncl ear. Note t h at the 
results of iKlev fc Dirksenl d2006l) and iFune et all (12014 ) 
may su pport the wide gap formation of IVarniere et al.l 
l| 20041) . IKlev fe Dirksenl d2006h also showed that the disc 
rotation has some eccentricity when the gap is extended 
to the 1:2 Lindblad resonance. T he eccentric gaps are 
formed for a large K (> 10 4 ) (e.g., IKlev fe Dirksenl 120061 : 
iFung et al .1120141 ). Such wide gaps by massive giant planets 
are beyond the scope of our one-dimensional disc model 
adopting the local approximation. 

In the above, we found that a larger Xd is required 
for a larger K (i.e., a massive planet) in order to repro¬ 
duce the minimum surface densities derived by DM13 and 
IVarniere et al.1 d2004l ). It should be noted that the propa¬ 
gation distance is not proportional to the parameter Xd- 
The propagating distance of waves is defined by the dis¬ 
tance from the wave excitation site to the angular momen¬ 
tum deposition site (i.e., *a). Since the one-sided toruqe is 
radially distributed (see Fig [Dr), the wave excitation site 
can be approximately given by the median of the distribu¬ 
tion, i.e., the point within the half of the one-sided torque 
arises. Such a excitation site shifts away from the planet with 


an increase of K (see Fig. II21FI . ID’Angelo fe Lubowl d20ld ) 
also showed this tendency that the peak of the excitation 
torque density shifts away from the planet as a deep gap 
is formed (Fig.15 in that paper), using hydrodynamic sim¬ 
ulations. Hence, because of this shift of the excitation site, 
th e propagating distance does not increase much as Xd with 
K. Goodman fc Rafikovl d200ll) showed that the propagating 
distance of waves decreases with an increase of the planet 
mass due to the non-linear wave damping. Hence, further 
studies are needed in order to confi rm whether this results 
given by Goodman fe Rafikovl d200ll) conflicts with ours be¬ 
cause of the shift of the excitation site. Furthermore, the 
non-linear wave damping would be weakened by the steep 
surface density gradi e nt at the gap edge, as pointed out by 
IPetrovich fe Rafikovl (|2012l ). In order to fix the parameter 
Xd, such a wave damping effect in the gap should be taken 
into account in future work. 


6 SUMMARY AND DISCUSSION 

We re-examined the gap formation in viscous one¬ 
dimensional discs with a new formulation. In our formula¬ 
tion, we took into account the deviation from Keplerian disc 
rotation and included the Rayleigh stable condition, consis¬ 
tently. We also examined the effect of wave propagation. Our 
results are summarized as follows. 

(i) The deviation from the Keplerian disc rotation makes 
the gap shallow. This is because of the enhancement of the 
shear motion and the viscous angular momentum transfer 
at the gap edges (see Fig. 0. 

(ii) For deep gaps, the deviation from the Keplerian disc 
rotation is so large that the Rayleigh stable condition is vio¬ 
lated. An enhanced viscosity dissolves such unstable rotation 
and makes it marginally stable (see Fig. 0. This effect also 
makes the gap shallower (see Fig. 0. 

(iii) To include the effect of wave propagation, we adopted 
a simple model where the position of the angular momen¬ 
tum deposition is parameterized by Xd- A large Xd indicates 
a long propagation length. The effect of wave propagation 
makes the gap wider and shallower (Fig. [9j. In a wide gap, 
the waves are mainly excited at the flat bottom, which re¬ 
duces the one-sided torque and the gap depth. For a suffi¬ 
ciently large Xd, the gap depth of our exact solution agrees 
well with the wide-limit gap and with the results of hydrody¬ 
namic simulations. At K = 1000, our model requires Xd ^ 6 
for the agreement (Fig. 1121) . In the case of instantaneous 
wave damping, on the other hand, our exact solution gives 
much deeper gaps than those of hydrostatic simulations. 

(iv) To check the validity of the large Xd, the gap width 
of our exact solution is compared with results of hydrody¬ 
namic simulations. For K = 1000, our exact solution with 
Xd ^ 6 has a gap width of 12 h p , whi ch is larger than thos e 
of DM13 (~ 8/ip). The gap widths of lVarniere et al] ((20041 ). 

4 In Fig. 1121 the exact solution have transition points from 
the realistic wave propagation (solid lines) to the unrealistic one 
(dashed lines) for each Xd . At the transition point of K, the ex¬ 
citation site defined by the median is equation to Xd . Fig. 1121 
shows that the excitation site moves away from the planet with 
an increasing K since the transition point of K increases with Xd • 
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on the other hand, are almost consistent with our exact so¬ 
lutions. Because of this uncertainty in the gap width of hy¬ 
drodynamic simulations, it is difficult to fix the preferable 
Xd by this comparison. 

(v) When the Rayleigh condition is taken into account, 
the deviation from the Keplerian rotation should also be 
included in order to keep the angular momentum conserva¬ 
tion. The Keplerian solutions with the Rayleigh condition 
give much shallower gaps, as shown in figures ISl and 1121 


In the case of instantaneous wave damping, the angular 
momentum deposition rate is given by A ex given by equa¬ 
tion (ESI- At far from a planet, the Rayleigh condition is al¬ 
ways satisfied because the planetary gravity is weak. Hence, 
the gradient of surface density is given by equation (1361) and 
the second derivative of surface density is given by 


d 2 Ins _ 4 KC 

dx 2 3* 5 


(Al) 


In future works, we need to det ermine the pref erable 


value of Xd■ Previous studies (e.g.. Takeuchi et al.l 

1991 

Korvcanskv & PaDaloizou 1996i: Goodman & Rafikov 

,2001; 

Dong et al. 201111 have investigated the wave propagation 


with no gap. As pointed out bv lPetrovich fe Rafikovl (20121 ). 
however, the gap structure can affect the wave damping. 
Since our result shows that the wave damping significantly 
affects both the gap depth and width, the wave damping 
should be treated accurately in both one-dimensional mod¬ 
els and hydrodynamic simulations for gap formation. 

Our simple model does not include the effect of the de- 
viation from Ke plerian disc rotation on the wave excitation. 
IPetrovich fc Rafikovl : 1* 2012 1 showed that a steep surface den¬ 
sity gradient modifies the excitation torque. Such a effect on 
the wave excitation should be included in future studies on 
the gap formation. Nevertheless, it is also considered that 
when the waves are mainly excited at the flat-bottom, such 
as for the wide-limit gap, the deviation of the disc rotation 
would not affect the wave excitation significantly. 

We also neglect the non-linearity of wave excitation, 
whereas the non-linearity cannot be neglected for large p lan- 
ets as A/p/M, > ( h p /R p ) 3 . According to iMivoshi et alJ 
(l999l) . the non-linearity makes the excitation torque small 
compared to the value for linear theory. This possibly leads 
to an additional shallowing effect. However, this effect would 
not significantly influence the gap depth since s m i n is scaled 
by only K in DM13’s relation (eq. [42]). 

The Rossby wave instability may be essential for the gap 
formation. In present study, we included only the Rayleigh 
condition. A more detail investigation including both the 
Rayleigh condition and the Rossby wave instability should 
be done in future works. 
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APPENDIX A: GAP MODEL BY 
TANIGAWA&IKOMA(2007) 

Following iTanigawa fe Ikomal (120071 4. we describe surface 
density structures of gap in a disc with the Keplerian ro¬ 
tation. In this section, we consider the gap structure in the 
case of instantaneous wave damping. The gap structure in 
the case with wave propagation is discussed in Appendix [B] 


The stability of the Rayleigh condition is checked by equa¬ 
tion (AH). Namely, when d 2 In s/dx 2 < — 1 in equation (AH, 
the surface density is described by the marginal Rayleigh 
stable state, instead of equation EH. We should point out 
that the second derivative given by equation m is used to 
determine the stability of the Rayleigh condition and does 
not affect the surface density distribution. Because of this, 
the surface density gradient is steeper and the shallowing 
effect of the Rayleigh condition is much higher than in our 
model. Using equation (30|l. we obtain the outer edge of the 
marginal Rayleigh stable region, x = x m , as 

4 \ 1/5 

3 CKj . (A2) 

In consideration of the continuity of the surface density dis¬ 
tribution, the surface density in the marginal Rayleigh stable 
region (x < x m ) is given by 

Ins = -|*m + \%m\x\ - ^x 2 , (A3) 

= -0.854K 2/S + 1.266K 1/f >| - 0.5x 2 , (A4) 

and the surface density for x > x m is given by equation EH- 
Since there is no torque density for x > 1.3 in this case, the 
minimum surface density s m i n is given by s(x = 1.3). Thus, 
we give s m in as 

Smin = exp (-0.854K 275 + 1.645K 1/S - 0.845) . (A5) 

Note that if x m < 1.3, the whole region of the gap is the 
Keplerian rotating part and s m i n is given by equation EH 
with x = 1.3. It should be also noticed that the angular 
momentum conservation is not satisfied at the flat bottom 
region with s = s m in of equation (IA5I) . as explained in the 
subsection 14. 3 I f and see also, Appendix G2>. 



APPENDIX B: GAP MODEL IN A KEPLERIAN 
ROTATING DISC WITH WAVE PROPAGATION 


Here, by assuming the Keplerian disc rotation, we derive 
gap solutions in the case with wave propagation. Following 
ITanigawa fe Ikomal (I20071 ). we consider the marginal condi¬ 
tion when the Rayleigh condition is violated. The angular 
momentum deposition rate is given by equation (1321) . Ignor¬ 
ing deviation in dfl/dR form the Keplerian and differentiat¬ 
ing, we give 


ds 

dx 


T hp_ 
3 w d 

0 


for < 1*1 < ZdT 

otherwise, 


(Bl) 
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Integrating equation m, we obtain the surface density 
without the Rayleigh condition as 


s = 



[(*<! +3^) “*] 


for |x| > Xd + 

for - ^ < |*| < Xd + ^ 
for |*| < *d - 

(B2) 


Equation (1B2I) does not satisfy the Rayleigh condition, 
especially for a large K (or a large T). First, the Rayleigh 
condition is violated near \x\ = *d + Wd/2h p because ds/dx 
of equation m is not continueous there. In order to make 
ds/dx continueous, the marginally stable condition (eq. 1291 1 
is used instead of equation m in the region where Xd + 
w d — (T7i p )/(3tOd) < |*| < *d + uid/Zhp. In addition, the 
marginally stable condition should also be used near |x| = 
Xd — Wd/2/ip for a large T. From equation (IB 1 1) . we find that 
the Rayleigh condition is violated from |x| = Xd — Wd/2h p 
to Xd + (wd/2/i p )(l — 6 /T) + 1. These two Rayleigh unstable 
regions are merged for T > 3(wd/h p ). In such large -K cases, 
since the marginally stable condition is used in the whole 
region of the angular momentum deposition, the minimum 
surface density is given by 



which is independent of K. This unrealistic minimum sur¬ 
face density does not satisfy equation E3 which is origi¬ 
nated by the an gular momentum conservation. As indicated 
bv iFung et al.l (2014J), the planetary torque should balance 
with the viscous angular momentum flux outside the gap for 
the angular momentum conservation. In our formulation, 
two terms in the right-hand side of equation m balance 
with each other in the bottom region (the left-hand side 
is negligibly small). However, the minimum surface density 
given by equation (1B3I) independent of K breaks down such 
a balance. Hence, equation (1B3I) also violates the angular 
momentum conservation. 

For a small K, equation (IB2I) is approximately valid 
because the Rayleigh unstable region near \x\ = Xd+Wd/2h p 
does not significantly affect s m i n . Substituting equation jB2j 
into equation (El, we give T as 


T = 3 


9A 3 

~KC 


1 - 


KCxd 


9[xl-(w d /2h p )*Y 


+ 1 


(B4) 


Substituting equation (lB2l) into equation (fB4l) . we obtain 
s m in for a small K as 


•Smin 


CK ( KCxd 

9A 3 y ~ 9 [a* - (w d /2 h p ) 2 f 



• (B5) 


APPENDIX C: SOLUTIONS FOR THE 
LINEARIZED EQUATION OF GAPS 

Here, we consider solutions to the following differential equa¬ 
tion: 

§-% = <?(*), (Cl) 

where g(x) is an arbitrary odd function of x. A general so¬ 
lution of this equation is given by a combination of homoge¬ 
neous solutions, e ± '^ x , and a particular solution. We seek 


a solution which satisfies the boundary conditions of y = 0 
at x = ±oo. Since g(x ) is odd in this case, the solution is an 
even function of x. For simplicity, we consider the solution 
for x > 0. Because of the symmetry of the equation, the 
solution of x < 0 can be obtained by inverting the sign of x 
in the solution for x > 0. 

A particular solution y v {x ) of equation (1C 111 can be 
given by 


Vr{x) ~ V5 [ 




POO 

/ s(*') 

J X 


e-^'dx 


-\/3i 


r 

J X 


g(x')e^ x dx 


(C2) 


In the case with instantaneous wave damping (see eq. [44]), 
g(x) is given by 


g(x) = 



3A 3 


for |x| > A, 


otherwise. 


(C3) 


Substituting equation (IC3I) into equation (IC2II . the particu¬ 
lar solution for x > A is obtained as 


Vp(x) = 


V3C 

18 


_\/3 3 

' x + 2 


-V3: 


: Ei{V 3*) - e^ x Ei(-V 3*) | 
(C4) 


wh ere Ei(ax) den otes the exp onential integral function (e.g., 
lAbramowitz fc Stegunlll965h . For x < A, the particular so¬ 
lution is given by 


y p (x) = a+e 
where a+ and a_ are defined by 


V3x , p -V3x , C 
+ + 9A 3 ’ 


nXU (C5) 


a± 


C_ 

6 


TK(T^A)-1(^ + 1),^ 


(C6) 


Using this particular solution given by equations (ICl and 
(IC5l) . we can obtain the general solution of equation El- 
Since this particular solution vanishes at x —» oo, the coef¬ 
ficient of the homogeneous solution e v ^ x is zero. The coeffi¬ 
cient of e ~'^ x , B, is obtained as 


1 dy p 

V3 dx 


(C7) 
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